c
c     This program is a test driver for the unsymmetric skyline solver.
c     This new version does not check the skyline results vs.
c     results from LINPACK.
c     
      integer          ndim, nndim
      parameter        (  ndim = 2500 )
      parameter        ( nndim = ndim*ndim+1 )
      double precision V(nndim), b(ndim), x(ndim)
      double precision normV, normx
      double precision bsave(ndim)
      double precision resid, residn, eps, epsmch, dmach
      double precision ops, sum, time(6), total
      double precision small, dens
      integer          n, job, band, job_setup
      integer          pd(ndim), nr(ndim)
      integer          NumProc
      logical          full
      character*1      afull
      parameter        ( small  = 1.0d-6 )

 
      full = .FALSE.
      job = 0
      epsmch = dmach('e')
      eps = 2.d0 * epsmch


c ___ Read matrix characteristics ___

      print *, 'Enter n'
      read( 5, * ) n
      print *, 'Is A full? ( Y/y for yes ):'
      read( 5, 999 ) afull
999   format(a)
      if( afull .EQ. 'y'  .OR.  afull .EQ. 'Y' ) then
	  full = .TRUE.
      else
	  print *, 'Enter bandwidth bw in % with respect to n',
     &             'bw <= 0 or bw > 100 generates random band'
	  read(5, *) band
      end if


c ___ Generate matrix and right-hand side ___
         job_setup = 3
         call dsetup( n, V, pd, nr, normV, full, band, b, job_setup )


c    ___ Save right-hand side for later checking ___

	 do i = 1, n
	    bsave( i ) = b( i )
	 end do


c    ___ Number of floating point operations. ___

	 if ( full ) then
 	    ops = (2.0d0*float(n)**3)/3.0d0 + 2.0d0*float(n)**2
	 else
	    sum = 0
            do i = 1, n
	       sum = sum + nr(i)*nr(i)
	    end do
	    ops = 2.d0*float(sum) + float(pd(n))/2.d0 +
     &		  2.d0*float(pd(n))
	 end if


c    ___ Factor the matrix and record the time ___

         t1 = second()
 	 call dskyfa( V, n, pd, info )
	 time(1) = second() - t1



c    ___ Solve using the factorization and time it. ___

	 t1 = second()
 	 if (job .EQ. 0) call dskysl( V, n, pd, b, 'N', info )
 	 if (job .NE. 0) call dskysl( V, n, pd, b, 'T', info )
	 time(2) = second() - t1


c    ___ Total time of factorization and solve. ___

	 total = time(1) + time(2)


c    ___ Compute a residual to verify results ___

	 do i = 1, n
	    x( i )    = b( i )
         end do

         call dsetup( n, V, pd, nr, normV, full, band, b, job_setup )

	 do i = 1, n
	    b( i )  = -bsave( i )
         end do


	 if ( job .EQ. 0 ) then
	    do j = 1, n
	       jnr = j - nr(j) - 1
	       jpd = pd(j) - nr(j) - 1
	       do i = 1, nr(j)
	          b(j) = b(j) + V( pd(j-1) + i )*x( jnr + i )
	       end do
	       do i = 1, nr(j)+1
	          b(jnr+i) = b(jnr+i) + V( jpd + i )*x( j )
	       end do
	    end do
	
	 else
	    do j = 1, n
	       jnr = j - nr(j) - 1
	       jpd = pd(j) - nr(j) - 1
	       do i = 1, nr(j)
	          b( jnr+i ) = b( jnr+i ) + V( pd(j-1) + i )*x( j )
	       end do
	       do i = 1, nr(j)+1
	          b( j ) = b( j ) + V( jpd + i )*x( jnr+i )
	       end do
	    end do
	
	 end if


	 resid    = 0.d0
	 normx    = 0.d0
	 do i = 1, n
	    resid = max( resid, abs( b( i ) ) )
	    normx = max( normx, abs( x( i ) ) )
         end do
	 residn    = resid / ( n * normV * normx * eps )


	 write(6,900)
 	 write( 6, 901 ) residn, resid, x(1), x(n)
900      format( / '     norm. resid      resid     ',
     &             '         x(1)          x(n)')
901      format( 1p4d16.8 )

	 call skyenvir( 'P', NumProc )
	 if ( full ) then
	     write(6,902) n
     &                  , NumProc
	 else
	     dens = pd(n) / ( float(n)**2 )
	     write ( 6, 905 ) n,
     &                        NumProc, 
     &                        pd(n), n, dens
	 end if
902      format(// '   Times are reported for matrices of order ', i5
     &           , '  on', i2, ' processor(s)'
     &         )
905      format(// '   Times are reported for matrices of order ', i5,
     &             '  on', i2, ' processor(s)',
     &        / '       ( density ', i7,'/',i4,'**2  = ', f5.2,' )' )
         write(6,903)
903      format(/ 6x,'factor',5x,'solve',6x,'total',5x,'mflops',7x)

	 time(3) = total
	 time(4) = ops / (1.d6 * total)
	 write(6,904) (time(i), i=1,4)
904      format( 4(1pe11.3) )

	 end
